function con = bcfunc_SK(s0, s1, el0, el1, mu,target_taf,chase_ta0,xyz)
%ta0 = atan2(sin(s0(1)-el0.raan), cos(el0.i)*cos(s0(1)-el0.raan)) - el0.aop;
ta0 = chase_ta0;
%[] calculate the true anomaly of the spacecraft

rtp0 = elements2rtp_xyz(el0, ta0, mu,xzy);
%[] perform coordinate transform from classical orbital elements to r,theta,phi

con(1) = s0(1) - rtp0(1);
%[] Initial R distance of the space craft

con(2) = s0(2) - rtp0(2);
%[] Initial theta position of the space craft

con(3) = s0(3)-rtp0(3);
%[] Initial phi position of the spacecraft initial position of the spacecraft in the x-y plane is constant

con(4) = s0(4) - rtp0(4);
%[] velocity in the r direction is constraint to initial condition

con(5) = s0(5) - rtp0(5);
%[] velocity in the theta direction is constrainted to initial condnition

con(6) = s0(6) - rtp0(6);
%[] velocity in the phi direction is constrainted to initial condnition

con(7) = s0(14);
%[] lagrange states of the thrust acceleration state is zero at initial time. therefore free.

ta1 = atan2(sin(s1(3)-el1.raan), cos(el1.i)*cos(s1(3)-el1.raan)) - el1.aop;
%[] calculate the true anomaly at final state

rtp1 = elements2rtp(el1, ta1, mu);
%[] Calculate the r,theta,phi coordinate position of the spacecraft

con(8)  = s1(1) - rtp1(1);
%[] Final r component constraint

con(9) = s1(2) - rtp1(2);
%[] final theta position constraints

con(10) = s1(3)-target_taf;
%[] Final lagrange multiplier of phi is zero, meaning the final state of the phi is free.

con(11) = s1(4) - rtp1(4);
%[] final r direction velocity constraint

con(12) = s1(5) - rtp1(5);
%[] final theta direction velocity constraint

con(13) = s1(6) - rtp1(6);
%[] final phi direction constraint

con(14) = s1(14) - 1;
%[] final lagrange state of the thrust accelreation = 1.
end